%This problem simulates the Circular Restricted Three Body Problem and
%experiments with possiblem solutions at various energies.

%clear the workspace
%clear;
%output long
format long;

%--------------------------------------------------------------------------
%----------------------NOTES-----------------------------------------------
%--------------------------------------------------------------------------

%NAME: CRTBPplayGround.m

%Dependencies:
%This program depends on 'CRTBP.m', 'librationPoints.m', 'jacobiConst.m'



%USE: Type 'CRTBPplayGround' from the MatLab command line. The user simply
%sets how long they want the integration to run.

%Below are several USER DEFINED CONTROL PARAMETERS.  The user should
%set these (or leave them set to their default values) and save the
%program.  

%OUTPUT:
%The output is a graph of the trajectory with the associated zero velocity
%curve.

%PURPOSE: To get a quick sample of the bahavior of the CRTBP.  Also tests
%the functions from the tool kit with random data.




%--------------------------------------------------------------------------
%--------------------USER DEFINED CONTROL PARAMETERS-----------------------
%--------------------------------------------------------------------------


    
    %NOTE:  The only parameter that has to be set by the user is the
    %integration time.  This is here:
        tf=20;    %default is tf=20.




    %Gravational constant
        G=1;

        
    %Set the value of mu:    
        %Let the program choose you a random value of mu between 0 and 1/2; 
        %mu=(1/2)*rand(1,1)


        %OR: Specify the value of mu you are interested in.  If you do this 
        %comment out the above.
        %The default is the earth/moon value.

        %for the earth moon system the value of mu is approx 0.012277471
        %mu=0.012277471;
        
        %ALSO:  if you want a random mu, but you want to skew toward the
        %likelihood of small mu you could use something like:
        mu=(1/2)*rand(1,1)
        %Or any power you like...
        
        
    %Set a random initial condition.  This is a vector of normally 
    %distributed numbers (so positive and negatives occur):   
    IC=(4/11)*randn(1,6)
    
    %Compute the Jacobi energy:
    C=jacobiConst(IC(1:3)', IC(4:6)', mu)
    
    
    %Compute the location of the five libration points as points in three
    %space.  These can be used to set the value of the Jacobi Constant.
    [L1, L2, L3, L4, L5]=librationPoints(mu);

    %Output location to the screen;
    outL1=L1
    outL2=L2
    outL3=L3
    outL4=L4
    outL5=L5
    %Now compute the Jacobi Energy at each collinear libration point

        %The velocity at any equilibrium is zero.
        v_equil=[0; 0; 0];

    %Compute the energy at each libration point
    CL3=jacobiConst(L3, v_equil, mu)
    CL1=jacobiConst(L1, v_equil, mu)
    CL2=jacobiConst(L2, v_equil, mu)
    CL4=jacobiConst(L4, v_equil, mu)
    %CL5=CL4 so no extra computation is needed
    
    

    

  



%------------------------SIMULATION----------------------------------------
%integrate the forward problem
    
    tspan=[0 tf];
    %numerical integration
    options=odeset('RelTol',1e-13,'AbsTol',1e-22);     %set tolerences
    [t,x]=ode113('CRTBP',tspan, IC,options,[],G,mu);



%------------------------------PLOTTING------------------------------------


%--------2d plotting-------------------------------------------------------

figure
hold on


%Contour Plot of the Energy Surface
   [X,Y]=meshgrid(-1.7:0.05:1.7);
    r1=((X + mu).^2+Y.^2).^(1/2);
    r2=((X + mu-1).^2+Y.^2).^(1/2);
    Z=2*[(1/2)*(X.^2 + Y.^2)+(1-mu)./r1+mu./r2];
    contour(X,Y,Z,[C,C],'r')


%plot the libration points
plot(L3(1), 0, 'ko', L1(1), 0, 'ko', L2(1), 0,'ko')
plot(L4(1), L4(2),'ko')
plot(L5(1), L5(2),'ko')

%plot the planets (primaries)
plot(-mu,0, 'r*', 1-mu, 0, 'g*')

%plot the forward trajectory
plot(x(:,1), x(:,2), 'b')


%-------3D PLOTTING--------------------------------------------------------

figure
hold on

%plot the libration points
plot3(L3(1), 0, 0, 'ko', L1(1), 0, 0, 'ko', L2(1), 0, 0, 'ko')
plot3(L4(1), L4(2), 0,'ko',  L5(1), L5(2), 0,'ko')

%plot the planets (primaries)
plot3(-mu,0, 0, 'r*', 1-mu, 0, 0,'g*')

%plot the forward trajectory
plot3(x(:,1), x(:,2), x(:,3), 'b')



